Analytical solution of the geodesic equation in Kerr-(anti) de Sitter space-times 
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The complete analytical solutions of the geodesic equations in Kerr-de Sitter and Kerr-anti-de 
Sitter space-times are presented. They are expressed in terms of Weierstrass elliptic p, (, and a 
functions as well as hyperelliptic Kleinian a functions restricted to the one-dimensional ^-divisor. 
We analyse the dependency of timelike geodesies on the parameters of the space-time metric and 
the test-particle and compare the results with the situation in Kerr space-time with vanishing 
cosmological constant. Furthermore, we systematically can find all last stable spherical and circular 
orbits and derive the expressions of the deflection angle of flyby orbits, the orbital frequencies of 
bound orbits, the periastron shift, and the Lense-Thirring effect. 



PACS numbers: 

I. INTRODUCTION AND MOTIVATION 

All observations in the gravitational domain can be ex- 
plained by means of Einstein's General Relativity. While 
for small scale gravitational effects (e.g. in the solar sys- 
tem) the standard Einstein field equations are sufficient, 
a consistent description of large scale obervations like the 
accelerated expansion of the universe can be achieved by 
the introduction of a cosmological term into the Einstein 
field equation 



kT, 



(1) 



where A is the cosmological constant which has a value 
of |A| < lfr 46 km~ 2 [J. 

Despite the smallness of the cosmological constant the 
question whether there might be measureable effects on 
solar system scales has attracted some attention. Within 
approximation schemes it was shown that all effects on 
this scale are too small to be detectable at present 0-0] ■ 
Nevertheless, there has been some discussion on whether 
the Pioneer anomaly, the unexplained acceleration of the 
Pioneer 10 and 11 spacecraft toward the inner solar sys- 
tem of apioncer = (8.47±1.33) x lfT 10 m/s 2 @, which is of 
the order of cH where H is the Hubble constant, may be 
related to the cosmological expansion and, thus, to the 
cosmological constant. The same order of acceleration 
is present also in the galactic rotation curves which as- 
tonishingly successfully can be modeled using a modified 
Newtonian dynamics involving an acceleration parame- 
ter omond which again is of the order of 10~ 10 m/s 2 @. 
Because of this mysterious coincidence of characteristic 
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accelerations appearing at different scales and due to the 
fact that all these phenomena appear in a weak gravity or 
weak acceleration regime, it might be not clear whether 
current approximation schemes hold. Therefore, it is dc- 
sireable to obtain analytical solutions of the equations of 
motion for a definite answer to these questions. 

There has been also some discussion if the cosmolog- 
ical constant has a measureable effect on the physics of 
binary systems, which play an important role in test- 
ing General Relativity. Although such an effect would 
be very small, it could influence the creation of gravi- 
tational waves @, §■ In particular, the observation of 
gravitational waves originating from extreme mass ratio 
inspirals (EMRIs) is a main goal of the Laser Interfer- 
ometer Space Antenna (LISA). The calculation of such 
gravitational waves benefits from analytical solutions of 
geodesic equations not only by improved accuracy, which 
is, in principle, arbitrary high, but also by the prospect 
of developing fast semi-analytically computation meth- 
ods @. Also, analytical solutions offer a systematic ap- 
proach to determine the last stable spherical and circular 
orbits, which are starting points for inspirals and, thus, 
important for the calculation of gravitational wave tem- 
plates. 

Finally, for a thorough understanding of the physical 
properties of solutions of the gravitational field equations 
it is essential to study the orbits of test-particles and light 
rays in these space-times. On the one hand, this is impor- 
tant from an observational point of view, since only mat- 
ter and light are observed and, thus, can give insight into 
the physics of a given gravitational field. On the other 
hand, this study is also important from a fundamental 
point of view, since the motion of matter and light can be 
used to classify a given space-time, to decode its structure 
and to highlight its characteristics. Furthermore, analyt- 
ical solutions give a possibility to systematically study 
limiting cases like post-Newton, post-Schwarzschild, or 
post-Kcrr expansions of geodesies and observables, which 
is also needed for a clear interpretation of the space-time. 
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Analytical solutions are especially useful for the anal- 
ysis of the properties of a space-time not only from an 
academic point of view. In fact, they offer a frame for 
tests of the accuracy and reliability of numerical inte- 
grations due to their, in pinciple, unlimited accuracy. 
In addition, they can be used to sytematically calcu- 
late all observables in the given space-time with the very 
high accuracy needed for the understanding of some ob- 
servations. In 1931 Hagihara [l(| first analytically in- 
tegrated the geodesic equation of test-particle motion 
in a Schwarzschild gravitational field. This solution is 
given in terms of the elliptic Weierstrass p function. 
The geodesic equations in Reissner-Nordstrom, Kerr, and 
Kerr-Newman space-times have the same mathemati- 
cal structure [ljj and can be solved analogously. For 
bound orbits in a Kerr space-time this has been elabo- 
rated recently [HI, [li|. The equations of geodesic mo- 
tion in space-times with non- vanishing cosmological con- 
stant exhibit a more complicated structure. Recently 
two of us found the complete analytical solution of the 
geodesic equation in Schwarzschild-(anti-)dc Sitter space- 
times based on the inversion problem of hypcrclliptic 
integrals [l4|, The equations of motion could be 

explicitly solved by restricting the problem to the 8- 
divisor, an approach which was suggested by Enolskii, 
Proninc, and Richter who applied this method to the 
problem of the double pendulum [ltjj . The mathemat- 
ical tool developed in these papers was also applied to 
geodesic motion in higher dimensional spherically sym- 
metric and static space-times [TtJ as well as to NUT-de 
Sitter and Plebahski-Dcmiahski space-times without ac- 
celeration [lj|. 

In this paper we extend the approach developed in 
[Til to the case of the stationary and axially sym- 
metric Kerr-(anti-)de Sitter space-times, thus generaliz- 
ing the results of both [ll| and [li|. We start with the 
derivation of the equation of motion for each coordinate 
dependent on proper time and decouple the equations for 
the r and 8 motion following an idea of Mino [lj| . Then 
we discuss possible types of test-particle orbits with a 
focus on the influence of the cosmological constant A. 
In section IIVI we explicitly solve the equations derived 
before and present for the first time the complete ana- 
lytical solution of the geodesic equation in Kerr-de Sitter 
space-time. After showing some chosen geodesies we de- 
rive the expressions of observables of particle and light 
trajectories. For bound orbits, the periastron advance 
and the Lcnsc-Thirring effect are given in terms of the 
fundamental orbital frequencies. 



where dr 2 = g^ v dx^dx u is the proper time along the 
geodesies and 



{&}= 2-9^ ( 5 p.9- + cUp„ - d„g pa ) (3) 

the Christoffcl symbol, in a space-time given by the met- 
ric 



dr 2 = (dt -asm 2 dd(f>) 2 - ^-dr 2 

X P A r 



A(> ° (adt - (r 2 + a 2 )dcpf - ^-d8 2 , (4) 
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where 



1 _ | r 2 ) ( r 2 + aj 



Ag = 1 + — COS 2 I 

a 2 A 

x = i + — , 

p 2 = r 2 + a 2 cos 2 



(5) 

(6) 

(7) 

(8) 



(in units where c = G = 1). This Boyer-Lindquist form of 
the Kerr-(anti)-dc Sitter metric describes an axially sym- 
metric and stationary vacuum solution of the Einstein 
equation and is characterized by M = 2M related to the 
mass M of the gravitating body, the angular momentum 
per mass a = J/M, and the cosmological constant A . 
Note that this metric has coordinate singularities on the 
axes 8 — 0, n and on the horizons A r = 0. The only real 
singularity is located at p 2 = 0, i.e. at simultaneously 
r = and 8 = J assuming a ^ 0. 

Analogously to the situation in Kerr space-time, we 
classify this form of the metric according to the number of 
(disconnected) regions where A r > 0, which depends on 
the parameters M, a and A. We speak of 'slow' Kerr-de 
Sitter if there are two regions and of 'fast' Kerr-de Sitter 
if there is one region where A r > 0. The limiting case 
where two regions are connected by a zero A r is called 
'extreme' Kerr-de Sitter. Other cases are not possible, 
what can be seen by a comparison of coefficients in A r = 
_A r 4 + (!_ ^a 2 )r 2 -Mr + a 2 = -^Ut=i( r ~^) whcrc 
r\j denote the zeros of A r . Fig. [1] shows the modification 
of regions of slow, fast, and extreme Kerr-de Sitter with 
varying A. 

We can identify four constants of motion, two corre- 
sponding to the energy per unit mass E and the angular 
momentum per unit mass in z direction L z given by the 
generalized momenta pt and p,p 



II. THE GEODESIC EQUATION 

We consider the geodesic equation 
d 2 x^ 







dr 2 



+ {£} 



dxP dx a 
dr dr ' 



(2) 



Pt = gut + 



gttj,t =: L z 



(9) 
(10) 



where the dot denotes a derivative with respect to the 
proper time r. In addition, a third constant of motion is 
given by the normalization condition 5 = g^x^x 11 with 
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FIG. 1: Regions of slow and fast Kerr-de Sitter (KdS) for 
different values of A. Left: the black curve corresponds to 
A = 0, the gray curve to A = — 10 -J . Below the line we 
have slow KdS and fast above. Right: A = 10~ 5 . The region 
bounded by the two curves corresponds to slow KdS, outside 
to fast. 



(5=1 for timelike and 5 = for lightlike geodesies. A 
fourth constant of motion can be obtained in the process 
of separation of the Hamilton- Jacobi equation 



2 dS_ _ ■ ■ dS dS 
dr ' 9 dx 1 dxi 



using the ansatz 
„ 1 



5t - Et + L z <j) + S r (r) + S 9 (9) 



(11) 



(12) 



If we insert this into (fTTj) we get 



5a 2 cos 2 9+A 6 



Sr 2 - A r 



dSr 
dr 



dSg 

89 

2 



Ag sin z 9 



(aE sin 2 9 - L z )' 



+ 7T ((r 2 + a 2 )E - aL z f , (13) 



where each side depends on r or 9 only. This means that 
each side is equal to a constant K, the famous Carter 
constant, [20| . 

From the separation ansatz (fT2"|) we derive the equa- 
tions of motion 



p 4 f 2 


= X 2 (EM 2 V) 2 - A r (5r 2 + K) , 


(14) 


p 4 9 2 


= Ag(K-a 2 S cos 2 9) -^ E ^\ 

sin 9 


(15) 


X 


a EM 

= —EM 2 ¥ 5— T , 

A r Ag sin 2 9 


(16) 


2 

X 


= — EM 2 V — T . 

A r Ag 


(17) 



where £M 2 P = (r 2 + a 2 )E - aL z and EMT = 
aE sin 2 9 — L z . The equations for r and 9 are coupled 
by p 2 = r 2 + a 2 cos 2 9. This difficulty can be overcome 
by introducing the Mino time A [l9[ which is related to 
the proper time r by ^jj = p 2 . For simplicity, we rescale 
the parameters appearing in eqs. (jl"4")) - (fr7|) such that they 



are dimcnsionlcss. Thus, we introduce 

r a - 1 . - , 

r=77. s =— , = 3 ' 



M 
L- 



M 
K 



M ' K M 2 ' 



(18) 



and accordingly 

A f = (1 - Af 2 ) (f 2 + a 2 ) - f , (A r = M 2 A f ) 
Ag = 1 + a 2 Acos 2 6>, (19) 
X = 1 + a 2 A. 

In addition, we can absorb (EM) in the definition of A 
by introducing 

^=§: «=J, 52 ^^' 7 = ™A. (20) 
Then the equations ()14 p - (fl7|) decouple and read 

2 



dr 
(W 

1 d(f> a 



— = R(f) = X V ~ Ar(S 2 f z + k) 



Q(9) =A e (K- S 2 a 2 cos 2 9) - 

sin 9 



P 



1 



X 2 d7 A f Ag sin 2 9 



1 dt 



f 2 + a 2 



Mx 2 c*7 
where, as before, 



P - — T. 

Ar Ag 



(21) 
(22) 
(23) 
(24) 



P = f 2 + a 2 ~ aD . 



(25) 
(26) 



T = a sin 2 0-V. 
In section [TVl we will explicitly solve these equations 

III. TYPES OF TIMELIKE GEODESICS 



Before solving the equations of motion derived in the 
previous section, we analyse the structure of possible or- 
bits dependent on the black hole parameters a, A and 
the particle parameters 5,E,L Z ,K. The major point in 
this analysis is that ([2T]) and (|2"2"|) imply R(f) > and 
Q(9) > as a necessary condition for the existence of 
a geodesic. Although we concentrate here on timelike 
geodesies light can be treated in the same manner and is 
in general easier to deal with. 

First of all, we state two short theorems about possible 
values of the Carter constant. The corresponding theo- 
rems for vanishing cosmological constant can be found in 

m- 

Theorem 1. If a geodesic lies entirely in the equatorial 
plane 9 = ^ or if it hits the ring singularity p 2 = then 
the modified Carter constant Q 
zero. 



K - X 2 (aE - L) 
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Proof. A geodesic lies entirely in the equatorial plane iff 
0( 7 ) = § for all 7. This implies that 9(6*) = (fp 2 = 
and with 



9 



K-X 2 (a-£>) 2 = QE~ 2 M- 



it follows Q = 0. If a geodesic hits the ring singularity, 
then there is a 7 such that f(7) = and 0(7) = -|. As 
R(f) > and 9(0) > for all r and 0, and in particular 
for f = 0, 6 = § , it follows 



j?(f = 0) = X 2 a 2 (a - Vf - a 2 K 



-« 2 Q 



and as above 



E 2 M 



i > 0. 



> => Q < 



□ 



Note that this theorem implies that Q = is a neces- 
sary condition for equatorial orbits, which is an impor- 
tant class found in many astrophysical objects like accre- 
tion discs and planetary systems. Note that the modified 
Carter constant Q depends on the cosmological constant, 
which also influences the next theorem. 

Theorem 2. For A > — all timelike and null geodesies 
have K > 0. In this case K = implies Q = and the 
geodesic lies entirely in the equatorial plane. 

Proof. A geodesic can only exist if there are values for 
f(7) and 0(7) with R(f) > and 9(0) > 0. From A > ^4 
it follows A g = 1 + a 2 A cos 2 > 1 - cos 2 > 0. If K < 
then (k — S 2 d 2 cos 2 0) < and 



9(0) = Ag(n — S 2 d 2 cos 2 



-T 2 < 



sin 2 



for all values of 0. Assume now K = 0. Consequently 



-T < 



9(0) = -6 2 d 2 cos 2 0A fl 



and 9(0) = only if cos 2 = and additionally T 



a sin 



V = a - V = 0. 



□ 



Since from observation the cosmological constant has 
a small positive value, the condition A > is always 
fulfilled. 

From these two theorems it is obvious that, while K 
originates from the separation procedure, the modified 
Carter constant Q has a geometric interpretation since it 
is related to possible values of the orbits. This relation 
will become more explicit in the following subsections. 

In the remainder of the section we will study the con- 
sequences of the two conditions R(f) > and 9(0) > 0. 



A. Types of latitudinal motion 

Geodesies can take an angle if and only if 9(0) > 0. 
Thus, wc want to determine which values of a. A, E, 



K, and E [0,7r] result in positive 9(0). For simplicity, 
we substitute v = cos 2 giving 



9(» = (1 + a A Kv){n - 5 2 a 2 v) 



X of {I - v) - 2aD 



V 2 



1 - v 



(27) 



Assume now that for a given set of parameters there ex- 
ists a certain number of zeros of 9(^) in [0,1]. If we 
vary the parameters, the position of zeros varies and the 
number of real zeros in [0, 1] can change only if (i) a zero 
crosses or 1 or (ii) two zeros merge. Let us consider 
case (i). is a zero iff 



9(i/ = 0) 



X 2 (a-V) 2 =0 



or 



aE± 



(28) 



(29) 



As v = 1 is in general a pole of <d(v) it is a necessary con- 
dition for 1 being a zero of 9(^) that this pole becomes 
a removable singularity. From (f2"7| it follows that this is 
the case for V = or, equivalently, L z = 0. Under this 
assumption we obtain 

Q{v = 1) = (1 + 5 2 A)(k- 5 2 d 2 ) forZ z = 0. (30) 

If we additionaly assume that A > (as in Thm. [2]) we 
can conclude that 9(^ = 1) = iff L z = and n = S 2 a 2 . 
Summarized, L z = aE ± and simultaneously L z = 

and k = 5 2 a 2 (assuming A > -^) give us boundary cases 
of the motion. 

Now let us consider case (ii) . If we exclude the coordi- 
nate singularities = 0, n or v = 1 the zeros of 9(z^) are 
given by the zeros of 

9„ = (1 - z/)(l + a 2 kv){n - 5 2 a 2 v) - \ 2 (a - V - av) 2 , 

(31) 

which is in general a polynomial of degree 3. Then two 
zeros coincide at x £ [0, 1) iff 



9„ = (1/ — x) 2 {a\v + a ) 



(32) 



for some real constants ai,ao- By a comparison of coef- 
ficients we can solve this equation for L z (x) and E 2 {x) 
dependent the remaining parameters a, A, and K. This 
parametric representation of values of L z and E 2 again 
correspond to boundary cases of the motion. Let us ad- 
ditionally consider the conditions for v — 1 being a dou- 
ble zero. With L z = and k = 5 2 a 2 (assuming A > ^) 
it follows 



dp 



a 2 x(x - 82) , 



(33) 



which is zero for x — $2 or, equivalently, E 2 



FIG. 2: Typical regions of different types of ^-motion. Here 
M = 2, A = 10~ 5 , and K = 3. The transition from slow 
Kerr-de Sitter (left, a = 0.4) to fast Kerr-de Sitter (right, 
a — 0.8) is continuous. 



For given parameters of the black hole a and A, we can 
use these informations to analyse the 9 motion of all pos- 
sible geodesies in this space-time. As a typical example 
for timelike geodesies consider Fig. [2] where the curves 
divide the half plane into four regions (a)-(d) which cor- 
respond to different arrangement of zeros in [0,1]. A 
geodesic motion is only possible in regions (b) and (d) be- 
cause in all other regions 0„ is negative for all v £ [0, 1]. 
Note that for the special case of k = 8 2 a 2 (assuming 
A > -^r) strictly speaking regions (b) and (d) are divided 
by L z = 0. However, in each region we have for L z > 
and L z < the same number of zeros in [0, 1] and, thus, 
the same type of motion. (More precisely, near L z = 
a zero v Q < 1 of Q(v) approaches 1, but does not cross 
it.) Therefore, in each region we put the parts above 
and below L z = together. The arrangement of zeros in 
the two regions (b) and (d) correspond to the following 
different types of motion in 9 direction (cp. Fig. [3]) 

• Region (b): 0„ has one real zero v max in [0, 1] with 
0„ > for v £ [0, z^max], i-e. 9 oscillates around the 
equatorial plane 6 = ^. 



Region (d): 0„ has two real zeros v m i n , z-Wx in [0, 1] 
with 0„ > for v £ [t'min, z-Wx], i-C 9 oscillates 
between arccos(±^/i/ m i n ) and arccos(±^/^ max ). 



The boundaries of region (b) arc given by L z = aE ± 



IK 

X 

and, therefore, the regions gets larger if grows, i.e. if 

K grows or a 2 A gets smaller. A change of a in addition 
causes region (b) to shift up or down. The dependence 
of region (d) on the parameters K, a and A is much more 
involved. The upper boundary of (d) is also the lower 
boundary of region (b) . The lower boundary is given in 
a complicated parametric form which makes it apparently 
impossible to determine an explicit connection between 
the form of region (d) and the parameters. However, the 
point where the upper and lower boundaries of region (d) 
touch each other is where is a double zero of 0„, which 



FIG. 3: Typical examples of Q(6) in the two regions (b) (left) 
and (d) (right). Allowed values of 6 are given by Q(6) > 0. 



is given by x = in E(x) and L z {x) from 



£7(0) 



la?_ 
2 



K - a 2 AK 



avK\ 



L z {0) 



1 a 2 - K X 



K X 



(34) 



The regions (b) and (d) are characterized in a simple way 
in terms of the modified Carter constant Q. As in region 
(d) 0^(0) < it follows that this region corresponds to 
Q < because of 0j,(O) = E 2 r i ■ In the same way we can 

conclude that region (b), where 0^(0) > 0, corresponds 
to Q > 0. 



B. Types of radial motion 

A geodesic can take a radial coordinate f if and only if 
R(f) > 0. The zeros of R are extremal values of f (7) and 
determine the type of geodesic. The polynomial R is in 
general of degree six in f and, therefore, has six possibly 
complex zeros of which the real zeros are of interest for 
the type of motion. As a Kerr-de Sitter space-time has no 
singularity in f = 0, 9 7^ f , we can also consider negative 
f as valid. However, r = is an allowed value of f(7) iff 



< R(Q) = x 2 (a 2 - aV) 2 - a 2 K 



Q 



E 2 r 



2 ©,(0). 



(35) 



It follows that f = can only be crossed iff Q < 0, which 
corresponds to region (d) of the 9 motion. In region (b) 
of the 9 motion where Q > a transition from positive 
to negative f is not possible. 

To clarify the discussion we introduce some types of 
orbits HH. 

• Flyby orbit: f starts from ±00, then approaches a 
periapsis f = r$ and back to ±00. 

• Bound orbit: f oscillates between to extremal n < 
f < T2 with —00 < n < ri < 00. 

• Transit orbit: f starts from ±00 and goes to =Foo 
crossing f = 0. 
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FIG. 4: Typical examples of R(r) for A > and for different 
numbers of real zeros. Allowed values of f are given by R(f) > 
0. 

All other types of orbits are exceptional and treated sepa- 
rately. They are either connected with the ring singular- 
ity p 2 — or with the appearence of multiple zeros in R, 
which simplifies the structure of the differential equation 
(PTj) considerably. Examples for the latter type are ho- 
moclinic orbits, cp. (22[. As large neg ative f correspond 
to negative mass of the black hole [23| , we will assign the 
attribute 'crossover' to flyby or bound orbits which pass 
from positive to negative f or vice versa. (By definition, 
a transit orbit is always a crossover orbit and, therefore, 
we will not explicitly state that.) Therefore, in region 
(d) of the 9 motion exists a crossover orbit, whereas all 
orbits located in region (b) of the 9 motion do not cross 
f = 0. 

For a given set of parameters we have a certain number 
of real zeros of R. If we vary the parameters this number 
can change only if two zeros merge to one. This happens 
at f = x iff 

R = (f- x) 2 {a A f A + a 3 f 3 + a 2 f 2 + a x f + a ) (36) 

for some real constants dj. By a comparison of coffi- 
cients we can solve the resulting 7 equations for E 2 {x) 
and L z (x) dependent on the remaining parameters a, A, 
and K. A typical result in slow Kerr-dc Sitter for small A 
including the results of the foregoing subsection is shown 
in Figs. [5] and [6l An analysis of the influence of each of 
the parameters a, A, and K is not done easily due to the 
complexity of the expressions for E 2 and L z . However, 
some typical examples for varying K are shown in Fig. [7] 
Examples of R(f) for different numbers of real zeros are 
given in Fig. [5] 

We discuss now the resulting types of orbits. For sim- 



plicity, we restrict ourselves here to the case of slow Kerr- 
dc Sitter although fast and extreme Kerr-de Sitter can 
be discussed analogously. This will be explicitly carried 
through in a future publication. 

For comparison, let us first study the situation in slow 
Kerr with A = 0. 

a. Case A = 0. We recognize five regions of different 
types of r motion. (Here we always assume r% < 

• Region (I): all zeros of R are complex and R(f) > 
for all r. Possible orbit types: transit orbit. 

• Region (II): R has two real zeros r 1; r 2 and R{r) > 
for f < r\ and r 2 < r. Possible orbit types: two 
flyby orbits, one to +oo and one to — oo. 

• Region (III): all four zeros n, 1 < i < 4, of R are 
real and R(f) > for r^h-i < f < r 2 k, k = 1,2. 
Possible orbit types: two different bound orbits. 

• Region (IV): again all four zeros of R are real but 
R(f) > for f < 7"i, r 2 < f < ?"3, and r A < f. 
Possible orbit types: two flyby orbits, one to each 
of ±oo and a bound orbit. 

• Region (V): R has two real zeros n, r 2 and R(f) > 
for n < f < T2- Possible orbit types: a bound 
orbit. 

Although there are the same number of real zeros the 
different orbit types in regions (III) / (IV) and (II) / (V) is 
due to the different behaviour of R when f — > ±oo. For 
A = the expression E 2 R = J2t=i a ^ ^ s a polynomial 
of degree 4 with 04 = E 2 — 1 which for f — > ±00 yields 
R(f) -> 00 if E 2 > 1 and R(f) -> -00 if E 2 < 1. 

Let us also analyse where we have crossover orbits. 
Region (II) is the only one which intersects region (b) as 
well as region (d) of 9 motion. As region (I) can only 
contain a transit orbit, which is by definition a crossover 
orbit, it can only intersect region (d). All other regions 
contain only region (b) of 9 motion and, therefore do not 
have any crossover orbits. The results of this paragraph 
together with the numbers of positive and negative zeros 
for each region are summarized in Tab. |TJ 

b. Case A > 0. Let us analyse now which regions 
change compared to the case A = 0. At first, we recognize 
that region (V) for A = merged with region (IV) and 
that region (III) becomes smaller for A > due to the 
shift of the separating E 2 = 1 line towards the left. A 
comparison of the possible orbit types for A > with the 
one for A = shows that in regions (I) and (II) there 
are no differences. However, these regions are slightly 
deformed (for small A) and a pair of parameters (E 2 , L z ) 
located in region (I) or (II) for A > may be located in a 
different region for A = 0. We consider regions (III) and 
(IV). Here again we assume < fj+i. 

• Region (III): all six zeros rj of R are real and R(f) > 
for f < r\,re < f and r 2 k < f < r 2 k+i for k 
1.2. Possible orbit types: two flyby orbits, one to 
each of ±00, and two different bound orbits. 
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0,0 0,5 1,0 



2,0 2,5 




FIG. 5: Regions of different types of f motion. Here a = 0.4 and K — 3. Left column A = 1CP 5 , middle A = and right column 
A = — 10™ 5 . Note that in the upper row for A > the vertical line has a maximum and for A < a minimum. For A = the 
vertical line is not related to a change of the number of real zeros but a switch from Iimf>->c» R(r) — oo to lmir->oo R(r) ~ — oo. 
Boundaries of the regions correspond to multiple zeros in R(r) and, therefore, to the occurrence of a stable or unstable spherical 
orbit. A detailed view including the results from the 8 motion can be seen in the lower row (note the rescaled axes). 



region 




+ 


range of r 


types of orbits 


Id 












transit 






lib 


1 


1 






2x flyby 






lid 


2 









flyby, 

crossover flyby 






Illb 





4 






2x bound 






IVb 


1 


3 






2x flyby, bound 






Vb 





2 






bound 








5 10 15 20 25 



TABLE I: Orbit types for A = 0. The + and - columns 
give the number of positive and negative real zeros of the 
polynomial R. Here the thick lines represent the range of 
orbits. Turning points are shown by thick dots. The small 
vertical line denotes f = 0. 



FIG. 6: Highlighted view on region (d) of 9 motion combined 
with r motion for a — 0.4, A = 10 , and K = 3. Regions 
(Id) and (lid) correspond to Q < and, thus, to crossover 
orbits whereas region (lib) corresponds to Q > 0. Note the 
differently scaled axes compared to Fig. [S] 



• Region (IV): R has four real zeros and R(f) > for 
f < r\,T2 < f < f3, f4 < f. Possible orbit types: 
two flyby orbits, one to each of ±oo and a bound 
orbit. 

Analogously to A = 0, regions (III) and (IV) only con- 
tain region (b) of 9 motion implying that there are no 
crossover orbits. Region (I) can only intersect region (d) 
because only transit orbits are possible. The remaining 
region (II) is the only one which intersects regions (b) 



and (d). 

We conclude that for E 2 > 1 the types of orbits are 
not noticeably changed, whereas for E 2 < 1 there are 
significant changes. In the former region (V) (for A = 0), 
which is now in region (IV), and in region (III) we have 
two additional flyby orbits which are not present for A = 
0. In a small vertical stripe left of E 2 = 1 there are even 
orbits which are bound for A = but reaching infinity 
for A > 0. In particular, it is independent of the value of 
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0,6 0,8 1,0 1,2 1,4 0,6 0,8 1,0 1,8 1,4 0,6 0,8 1,0 1,8 1,4 0,6 0,8 1,0 1,2 1,4 



FIG. 7: Deformation of regions of r motion for a = 0.4, A = 10~ 5 , and different K. From left to right K = 0, K = 1, K = 5, 
and K = 9. For K = region (d) of 9 motion vanishes and region (b) reduces to the line L z = aE, which, therefore, is the 
only allowed parameter region and corresponds to Q = and 6 = -| , see Thm. [2] Around this line L z — aE additional parts 
of regions (III) and (IV) appear not present for K > 0. 



region 
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range of r 


types of orbits 


Id 












transit 






lib 
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2x flyby 






lid 
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flyby, 

crossover flyby 






Illb 
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2x flyby, 2x bound 


IVb 


1 
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2x flyby, bound 







TABLE II: Orbit types for small A > 0. For the description 
of the +, — and range of f columns see Tab. [I] 
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range of r 


types of orbits 


Id 




1 






crossover bound 


1 






lib 










2x bound 


2 


2 






lid 




1 






bound 

crossover bound 


3 






Illb 











2x bound 


4 






IVb 










3x bound 


2 


4 







TABLE III: Orbit types for A < 0. For the description of the 
+ , — and range of f columns see Tab. [I] 



E if a geodesic may reach infinity as expected from the 
repulsive cosmological force related to A > 0. 

Note that for huge A the separation in regions (I) to 
(IV) is no longer possible because the repulsive cosmo- 
logical force becomes so strong that bound orbits are no 
longer possible. In this case we have only two regions, 
one with two real zeros corresponding to two flyby or- 
bits and one with only complex zeros corresponding to a 
transit orbit. 

All orbit types for small A > are summarized in 
Tab.HH 

c. Case A < 0. Here region (V) from A = merges 
with region (I) and region (III) becomes larger for A < 
due to the shift of the E 2 = 1 line to the right. Compared 
to the situation for A = the possible orbit types in 
region (III) do not change, but a set of parameters located 
there may be located in a different region for A = 0. Let 
us examine the remaining regions. 

• Region (I): R has two real zeros n < r-i and R(f) > 
for ri < f < r 2 . Possible orbit types: bound 
orbit. 

• Region (II) (and (III)): R has four real zeros with 
R(f) > for r 2 fc-i < f < r 2 k, k = 1,2. Possible 
orbit types: two different bound orbits. 

• Region (IV): all six zeros of R are real and R(f) > 



for r2fc-i < r < r2k, k = 1,2,3. Possible orbit 
types: three different bound orbits. 

Concerning crossover orbits regions (III) and (IV) again 
only contain region (b) of 6 motion. Also region (II) 
intersects both (b) and (d) whereas region (I) can only 
contain region (d) of 6 motion. 

Summarizing, the types of orbits significantly change 
if E 2 > 1. The transit orbit in region (I) for A = 
is transformed to a bound orbit for A < as well as the 
flyby orbits in regions (II) and (IV). Although region (V) 
for A = merges with region (I) for A < 0, the types of 
orbits do not change there. In general, because of R — > 
— oo if f — > ±oo we can not have orbits reaching f = ±oo 
at all as expected due to the attractive cosmological force 
related to A < 0. 

All orbit types for A < are summarized in Tab. IIIII 



IV. ANALYTIC SOLUTIONS OF THE 
EQUATIONS OF MOTION 

We will now analytically solve the geodesic equation in 
Kerr-de Sitter space-time (|2ip - (|24j) . Each equation will 
be treated separately. 
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9 motion 



We begin with the differential equation ([22 

■ d6 \ 2 ~2tt 



v 2 T 

9 = A e (K-5 2 a 2 cos 2 ^ x 



sin 2 (9 ' 



which can be simplified by the substitution v = cos 2 9 
yielding 



1 ( dv 



4 \d^i 



(37) 



where 0„ is defined as in (|31[) . This differential equation 
can be solved easily if vQ v has a zero with multiplicity 2 
or more. In this case (l37l) can be rewritten as 



4(7 - 70 ) 



(38) 



where 70 and v$ arc initial values, P 2 is a polynomial with 
maximum degree 2, and is a zero of vQ v with multiplic- 
ity 2j or 2j ' + 1, j = 1, 2. The integral on the right hand 
side can then be solved by elementary functions [24J. As 
in this case the explicit expression provides no further 
insight and some case distinctions would be necessary we 
skip the solution procedure. 

If vQ v has only simple zeros the differential equation 
(|37[) is of elliptic type and first kind and can be solved in 
terms of the Weierstrass elliptic function p. In contrast 
to the r motion considered in the next subsection, this 
structure does not simplify if we consider only light with 
6 = 0. To obtain a solution we transform vQ v to the 
Weierstrass form (4y 3 — g 2 y — 33) for some constants g 2 
and 33: First, we substitute v = giving 



(39) 



where 



B 6 := £ 3 (k - x 2 (d - Vf) + £ 2 (a 2 (/iA - S3) - k 

+ 2 X 2 a(a - V)) + a 2 (<5 2 (l - Aa 2 ) - X 2 - Ak)£ 
+ 5 2 d 4 A (40) 

=:£«*€*. 

i=l 

Second, we substitute £ = ^- (4y — ^) yielding 



where 



1 { dy\ 



O2 aifl3 
• 9 2 = 12 _ — 



(41) 



Q3 = — aiao&s dntto a% . (42) 

y 48 16 3 216 2 y ' 



1 .2 1 3 



The differential equation ([4T[) is elliptic of first kind, 
which can be solved by 

y{j) = p( 2 7 - 70,in; 92, 93) ■ (43) 
Accordingly, the solution of (|22|) is given by 



0(7) = arccos ( ± 



C63 



4p(27-7e,i„;ff 2 ,33) - f- 



(44) 



where 701 



270 + 11 



dy_ 



^iy' 3 -g 2 y'~g 3 



with j/o = 



4cos^(e ) T2 d e P en ds on the initial values 70 and 9q 
only. The sign of the square root depends on whether 
0(7) should be in (0, ^) (positive sign) or in (f ,7r) (neg- 
ative sign) and reflects the symmetry of the 9 motion 
with respect to the equatorial plane 9 = -|. If the mo- 
tion is located in region (b) from the previous section 
this implies that the two solutions have to be glued to- 
gether along 0(7) = \ if the whole 9 motion should be 
considered. 



B. 



r motion 



The differential equation that describes the dynamics 
of r (EU 



d r 



R = X 2 F 2 - Ar(d 2 f 2 



is more complicated because R is a polynomial of a degree 
up to 6. If R has a zero of multiplicity 4 or more, or if 
R has two zeros of multiplicity 2 or more, the differential 
equation (f2Tj) can be written as 



7-70 



df' 



r UU(r'-n) ji V^W) 



(45) 



where 70 and fo are initial values, P 2 is a polynomial 
with maximum degree 2, fj arc zeros of R with multi- 
plicity 2ji or 2ji + 1 where ji = 1,2, and k = 2 if there 
are two zeros of multiplicity 2 or more and k = 1 else. 
The integral on the right hand side can then be solved by 
elementary functions [24| . As the explicit expression pro- 
vides no further insight and some case distinctions would 
be necessary wc skip the solution procedure. 

If we consider null geodesies, i.e. S = 0, R is in general 
of degree 4 and the differential equation ([21) is of elliptic 
type and first kind. Then it can be handled using the 
method presented in the foregoing subsection: With the 
substitutions f = + Tr, where ffj is a zero of R, and 

£ = ~ (4^ 
at a form fll]). 



y), where bi 



we arrive 



This can then again be solved in terms 
of Weierstrass elliptic functions. The result is 



r( 7 ) = 



4p( 7 



7f,in;ff2,53) - ^ 



(46) 
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where 



7f 



= 7o 



f 
•> 1/1 



dy' 



with 



2/0 = 



b 3 



y/4y' 3 -g 2 , r y'-g 3 ,T 

4(fo-r R ) if depends only on the initial values 70 and 
fo and (72 j 53 are defined as in (|42[) with a j = 6, . 

The differential equation (|21|) is also of elliptic type 
but of third kind if R has a double or triple zero f 1 . In 
this case (|2~T1) reads 



7-7o 



f/r' 



(r-rx)v^(F) ' 



(47) 



where P4 is a polynomial of degree 4. This equation can 
be solved for f (7) analogous to the method which will be 
presented in subsection IIV CI 

If we consider particles, i.e. (5=1, and assume that R 
has only simple zeros the differential equation (f2~Tj) is of 
hyperelliptic type. It can be solved in terms of derivatives 
of the Kleinian a function with the method developed in 
14 1 . For this, we have to cast (|2ip in the standard form 
by a substitution f = ±— +fn with a zero r/j of i?. This 
yields 



du\" 
u— = c 5 R u , 

djj 



where 



(6-i)\ du^-V 



(48) 



(Tr) ■ (49) 



The sign in the substitution has be chosen such that the 
constant C5 is positive and, therefore, depends on the 
choice of rjj and the sign of A. The differential equation 
(|48|) is of first kind and can be solved by 



^2 I v / C^7-7ran , 



(50) 



where 7 f , in = + / 



00 u'du' 



and u = ± (r - r R ) 



depends only on the initial values 70 and f$. Here / is the 
function that describes the £>-divisor, i.e. a ((f(x), x) 1 ) = 
0, cp. The radial distance f is then given by 



f( 7 ) = T ^ f /(VCBT - 7f, in) 
CT 1 \ v/C57-7f,in , 



Hi: 



(51) 



where the sign depends on the sign chosen in the sub- 
stitution r = ± — + tr, i.e. is such that C5 in (j4"9")l is 
positive. 



motion 



This equation can be splitted in a part dependent only on 
f and in a part only dependent on 9. Integration yields 



30 = X 



70 A K7) 



G?7 



7 



, A e(7) sin 0(7) 
Td6' 



dPdf' 



AfVR Je A e > sin 2 6'y/& 



where we substituted r = r (7), i.e. 



^7 



(52) 



i?, in the first 



and 9 = 0(7), i.e. |^ = V0, in the second integral. 
We will solve now the two integrals in (|52p separately. 



The 9 dependent integral 



Let us consider the integral 



V] d0' 



On 



An' sin 



which can be transformed to the simpler form 



/e = T 2 



a — T> — av' 



A^(i-!/')vve7 



-.du' 



(53) 



(54) 



by the substitution v = cos 2 9, where Q v is defined in (|3"Tj) 
and A^ = 1+a Av. Here we have to pay special attention 
to the integration path. If 9 G (0, ^] we have cos 6* = 
+^/v but for 9 £ 7r) it is cos# = — Accordingly, 
we first have to split the integration path from 9q to 9 
such that every piece is fully contained in the interval 
(0,^] or [§,7r) and then to choose the appropiate sign 
of the square root of v. In the following we assume for 
simplicity that cos# = 

Analogous to subsection IIV Al the integral Ig can be 
solved by elementary functions if vQ v has at least a dou- 
ble zero [24[ ■ If v®v has only simple zeros Ig is of elliptic 
type and of third kind. If this is the case, the solution to 
Ig is given by 



«:5 



C(vi)(v - v ) + l0£ 



<t(v - Vi) 



2?t ik, 



• (a 3 A( x - aAV)(S a + S i2 ) + V{5 a + 8 U )) \ (55) 



We treat now the most complicated equation of motion 
in Kerr-de Sitter space-time, namely the equation for the 
azimuthal angle ([23]) 



d"f 



Ar 



T 



A, 



where the constants m are defined as in subsection IIV A( 



pM = ft 



1 T,2 



Aa 3 = p(«a), p(v3) = ft + f = 



u = v(j) = 2-f-^/g^n with 7 ejin as in (03]) and w = w(7o)- 
The integers fc, correspond to different branches of log. 
The details of the computation can be found in appendix 

m 
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2. The r dependent integral 
We solve now the first, f dependent integral in (|52p 

t? - aD) df' 



Ir : = 



a if' 2 + a z 



(56) 



Analogous to subsection IIVBI this integral can be solved 
by elementary functions if R has a zero with multiplicity 
4 or more or two zeros with multiplicity 2 or more [24j . 

If we consider light, i.e. 6 = 0, R is in general of 
degree 4 and I r is of elliptic type and third kind. In 
this case it can be solved analogously to Ig. The same 
substitutions f = | + r R and £ = (4j/ — -j) as in 
subsection II V Bl for the case J = 0, a subsequent partial 
fraction decomposition, and the final substitution y = 
p(v) result in 



dv 

, P(v) - yi 



d (f| + a 2 - aP) 

^f—f ft 



dv , 
(57) 



where j/j are the four zeros of A J/ ( ? ), 63 defined as in 
(|46|) . and are the coefficients of the partial fractions 
dependent on the parameters and tr . The four functions 
fi(v) — (p(v) — have simple poles in Vn,Vi% with 

p{vii) = 2/; = p( v i2) and have to be integrated wih the 
method presented in appendix [X] Then I r is given by 



4 2 
r _ \ "* \ ^ Cj 



log cr(w - Vij) 



C(vij)(v - v ) + log a(v - Vij) 
"I + a 2 - aX>) 



■(«-«o), (58) 



wo = w(7o) the result is 

■|ci(«;-t«o) + Co(/H-/(t«o)) 



E 



/C5KI 

/ Km 



1 .t(VK+H) 1 a(^+( Wo )) 

2 ° g a(W-(«;)) 2 ° g <7(W- («*,)) 



(/M - f(u> ),w- w 




(59) 



where the constants Cj are the coefficients of the par- 
tial fractions, Ui are the four zeros of Af—±i/ u+ f R , uq = 
±(r— r/j) — 1 , and C5 and R u as in (gSJ). The functions 

are defined by W ± (w) := (f(w),wY ~ 2 J * g?z, where 
the points uf = (ui, ± */ ) on the Riemann surface of 
y 2 = R u are the pole Uj located on the positive and the 
negative branch of the square root. 

D. t motion 

The equation for t (JM|) 



— = X 2 M 
cry 



r 2 + a 2 a 
— -P-—T 

A* Ag 



is as complicated as the equation for cj) motion. An inte- 
gration yields 





- p 


= X 2 M 




= X 2 M 


/' 

Jr 


= X 2 M 


Ir- 



Prf7 - 



(f 2 + a 2 )P 



A 



R 



dr 



-Td-y 



aT 



:d6 



^e 1 



(60) 



Because we already demonstrated the solution procedure, 
we only give here the results for the most general cases. 

If v<d u in (|37|) has only simple zeros the solution of the 
9 dependent part is given by 



where v = v(j) = 7 — jf.in, vq = u(7o) with jf.in as in 
(|46| . In the same way I r can be solved if R has a double 
or triple zero. 

If we consider particles, i.e. 6=1, and assume that R 
has only simple zeros, I r is of hypcrclliptic type and third 
kind. The details of the solution method can be found in 
appendix iBl but we give an outline here: First, we trans- 
form analogously to section HVB I to the standard form by 
f = ±1/ u+tr with a zero tr of R. Afterward we simplify 
the integrand by a partial fraction decomposition which 
allows us to express I r in terms of the canonical holo- 
morphic differentials dz (|B10|) and the canonical differ- 
ential of third kind dP(xi,x 2 ) (|B7|) . These differentials 
can then be expressed in dependence of the normalized 
Mino time 7. If we define w = w(j) = v/c57 — 7f,in and 



~ 2 A 

Ig = a 3 (v -V )-Y] * r [£(Vi)(v ~ V ) 

4p'(ui) 

+ log(r(v-Vi)-loga(v -Vi)] (61) 

where 03 is defined as in (|3D|l. p{v\) = ff — jd 2 Aa 3 = 
p(v 2 ), and v = v(j) = 2j - jg )in , v = v(j ) as in (JS3J). 

If we consider light, i.e. 6 = 0, the solution for the f 
dependent part is given by 



&3 j _ \ ^ \ ^ Cj 



C(vij){v - v ) + log a(v - Vij) 



loga(w - 1\ 



(fR + d 2 )(f 2 R + d 2 -dV) 



A r=rR 



(v - v ) , 

(62) 
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where 63 is defined as in (|46|) . C, are the coefficients of 
the partial fractions, p(vn) = yi — with the four 

zeros y l of A y(f) , y and v = v(~f) = 7- 7^, = u (7o) 
as in (gB]). 

If i? has only simple zeros and we consider timelike 
geodesies 8=1 the solution of the f dependent part is 
given by 





E 



2.i 



C 1 (w-wo) + C (f(w)-f(w )) 



1 oQ^+H) 1 , a(VF+( Wo )) 

2 ° g a(W-(u;)) 2 ° g (T(VF-(u;o)) 



(/M - f(w ),w- w ) 




(63) 



where the notation is as in (|59|) and Co,Ci,C2,i are the 
coefficients of the partial fractions. 



V. DISCUSSION OF SOME GEODESICS 

In Sec. IIIII wc discussed the general features of the dif- 
ferent types of timelike geodesic motion in Kerr-de Sitter 
and Kerr-anti-de Sitter space-time. With the analytical 
solution derived in section Hvl at hand we want to discuss 
now some chosen geodesies. 

We start with orbits which highlight the influence of 
A on the geodesies. From the results of section Mil we 
conclude that for A > there are four parameter regions 
where the changes compared to A = are most obvious. 
The first two are the regions (III) and (IV) with E 2 < 
1, where we have additional fly by orbits not present for 
A = 0. Third and fourth, the shift from region (V) of 
A = to region (II) of A > for E 2 = 1 - e, e > small, 
and the shift from region (III) of A = to region (IV) of 
A > 0, again for E 2 = 1 — e are most interesting as the 
(outer) bound orbit becomes a flyby orbit. A plot of the 
corresponding orbits can be found in Fig. [8] 

An important feature of geodesies in stationary ax- 
isymmetric space-times is motion of the nodes where the 
orbit of a test-particle or light intercepts the equatorial 
plane. This motion is caused by the got components of 
the space-time metric and known as the Lense-Thirring 
effect. In the weak field regime it becomes visible by a 
precession of the orbital plane, cp. Fig. [9] for an obvi- 
ous example. This orbital precession has been confirmed 
within an accuracy of about 10% by the LAGEOS (Laser 
Geodynamics Satellite) mission [2a| . j37| . 

Let us also discuss some exceptional orbits related to 
multiple zeros of R, i.e. spherical orbits with constant 
r and orbits asymptotically approaching a constant r. 
There are two types of spherical orbits: stable and un- 
stable. Stable spherical orbits with f (7) = f occur if 
radial coordinates adjacent to fo are not allowed due to 
R(f) < 0, which happens if fo is a maximum of R. Un- 
stable spherical orbits with ^(7) = fo are trajectories 



(a)E 2 = 0.9, L z = -1 
(region (IVb)) 



(b)E 2 = 0.94, L z = 0.6 
(region(IIIb)) 





= 0.97, L z = -1 
(region (lib)) 



(d)E 2 = 0.97, L z = 1 .. 
(region (IVb) 



FIG. 8: Flyby orbits for A = 1(T 5 , a = 0.4, and K = 3. 
For A = the parameters of (c) would have been located in 
region (Vb) and the parameters of (d) in region (Illb). Light 
grey cones correspond to extremal 9 and dark grey spheres to 
horizons. 




FIG. 9: Precession of the orbital plane for a bound orbit 
with A = 10" 5 , a = 0.4, K = 4.5, L z = -0.5, and E 2 = 0.96 
(region (Illb)). As in Fig. [8] the cones and spheres correspond 
to extremal 8 and horizons. 



where radial coordinates f in the neighbourhood of ?~o 
with f < fo or f > fo are allowed. Therefore, these or- 
bits are related to a minimum or to an inflection point of 
R. If fo is an inflection point, an asymptotic approach 
to ?~o is only possible from one side of fo whereas this is 
possible from both sides if fo is a minimum of R. Asymp- 
totic orbits can also be devided into two types: unbound 
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and bound. The latter case corresponds to orbits which 
approach for both t — ¥ oo and t — > — oo a spherical or- 
bit. Bound asymptotic orbits are also known as homo- 
clinic orbits. If the asymptotic orbit is unbound it reaches 
f = oo for either t — > oo or t — > — oo. 

For asymptotic bound or unbound trajectories corre- 
sponding to an unstable spherical orbits the equations of 
motion simplify considerably. In this case the equation 
for f(7) as well as the f dependent integrals in the 4> 
and t equations are of elliptic type and can be solved in 
terms of Weierstrass elliptic functions, see (|47f and (f58|) . 
Note that these solutions are not limited to the case of 
equatorial circular orbits but are valid for all types of 
asymptotic orbits and, thus, generalize the analytical so- 
lutions for homoclinic orbits in 22[ not only to Kerr-de 
Sitter space-time but also to arbitrary inclinations. 

From all spherical orbits the Last Stable Spherical Or- 
bit (LSSO), in particular, the Innermost Stable Circular 
Orbit (ISCO) in the equatorial plane are of importance as 
they represent the transition from stable orbits to those 
which fall through the event horizon. The corresponding 
multiple zero of R appears at the boundaries of the dif- 
ferent regions of r motion, cp. Fig. [5] Because necessarily 
K = x 2 E 2 (a — T>) 2 for equatorial orbits, from this we can 
determine the LSSO for given a, A, K and the ISCO for 
given a, A by solving first 




(a) a = 0.1, A = 1CT 



K : 



(b)a = 0.4, A = 10" 



K : 



FIG. 10: (a): Last Stable Spherical orbit at r(-y) = 
2.98982782315. The corresponding parameter values are (ap- 
proximately) E 2 = 0.88830740016 and L z = 0.00154804028 
(lower left corner on boundary of region (Illb)). (b): Stable 
spherical orbit at f(7) = 0.14 within the inner Cauchy hori- 
zon. Here E 2 = 124.76756414675 and L z = 5.879796033955 
(region (lib)). 




R(f) = 0, ^(r) = 0, and ^J(f) = (64) 



for f > f/j with the event horizon f^. The solutions are 
limiting cases of the LSSO or ISCO and are given by 
the corner points on the borders of region (III) of the f 
motion (as corners on other boundaries correspond to f < 
fh). From the results of (|64|) we search for the smallest 
possible double zero f which is a maximum. In the case 
of the ISCO in the equatorial plane we are now done. 
For the LSSO, we have to check in addition whether the 
corresponding values of E 2 (f) and L z (f) (given by ([36]) ) 
are located in an allowed region of the 9 motion. If this 
is the case, we found the LSSO. If not, we can determine 
the LSSO as the intersection point of the boundary of 
region (III) with a boundary of an allowed 8 region. Note 
that it is not possible to determine an LSSO (for given 
a, A, and K) if there is no spherical orbit at all outside 
the event horizon which happens if no boundary of the f 
motion is located in an allowed region of the 9 motion. 
As an example, this is the case for A = 10~ 5 , a = 0.1, 
and K = 0.1. Also, the LSSO is identical with the ISCO 
if it is given as an intersection point with the boundary 
of region (b) of the 9 motion. For examples of spherical 
orbits see Figs. [TU] and [TTJ Note that within the event 
horizon there may be additional stable spherical orbits. 



FIG. 11: Unstable spherical orbit with r(-y) = 1.5 (left) and 
asymptotic approach (right). Here a — 0.3, A — 10~ 5 , K — 
2, E = 0.9720311146 and L z = 2.5677957370 (boundary of 
region (Illb) ) . 



VI. ANALYTIC EXPRESSIONS FOR 
OBSERVABLES 

For the understanding of characteristic features of 
space-times by measurements of geodesies in that space- 
time it is crucial to identify certain theoretical quantities 
of observables. For flyby orbits, this can be the deflec- 
tion angle of the geodesic whereas for bound orbits it is 
of interest to determine the orbital frequencies as well as 
the periastron shift and the Lense-Thirring effect. 

Let us first consider flyby orbits. The deflection angle 
of such an orbit depends on the two values 7^ of the 
normalized Mino time for which ^(7^) = 00. These are 
given by 



1 



7 ± 

/ 00 



for the two branches of 
late the values of 8 and 
this flyby orbit by 8 ± = 



udu 



(65) 




Therefore, we can calcu- 
00 that are taken by 
and 0± = 0( 7 ± ). The 
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deflection angles arc then given by A9 = 9 + — 9 and 
Acj) = (j) + - cj)-. 

For bound orbits we can identify three orbital frequen- 
cies Q r , £lg and fi<£ associated with the coordinates r, 9 
and <j). The precessions of the orbital ellipse, which in 
the weak field regime can be identified with the perias- 
tron shift, and the orbital plane, which in the weak field 
regime can be identified with the Lense-Thirring effect, 
are induced by mismatches of these orbital frequencies. 
More precisely, the orbital ellipse precesses at f2^ — Sl r 
and the orbital plane at tt^ — Qg. 

Let us consider the orbital frequency f2 r . For bound 
orbits the coordinate f is contained in an intervall f p < 
f < f a with the peri- and apoapsis distances f p and f a . 
The orbital period u>f defined by f(7 + cjf) = f('y) is then 
given by a complete revolution from f p to r a and back 
(with reversed sign of the square root) to r p , 



dr 

Vr 



(66) 



The orbital frequency of the r motion with respect to 
7 is then given by — . For the calculation of O r , which 
represents the orbital frequency with respect to t, we need 
in addition the average rate T at which t accumulates 
with 7. This will be determined below. 

For the calculation of the orbital frequency fig we have 
to determine the orbital period ug such that 9(j + ujg) = 
0(7). The 9 motion is likewise bounded by 9 m [ n < 9 < 
#max for two real zeros 9 m - ln , 9 max £ (0, tt) of 8 and, there- 
fore, 



LJg 



d9 



(67) 



Again, the orbital frequency of the 9 motion with respect 
to 7 is given by 

The orbital periods of the remaining coordinates t and 
4> has to be treated somewhat differently because they 
depend on both f and 9. The solutions t(j) and 
consist of two different parts, one which represents the 
average rates T and at which t and <j> accumulate 
with 7 and one which represents oscillations around it 
with periods uj f and tug [TH, (2(| . The periods T and Y$ 
can be calculated by [HI 



2 

UJf 



f « (f 2 + a 2 )P 



Y* = — 



A r VR 
aP 



dr 



2 

U9 Jfl 



dr — 



dTd9 
AgVQ 
Td9 



(68) 
.(69) 



A r \fR ' Je mln A e sin 2 9V@ 

The orbital frequencies Q r , and are then given by 

(70) 



2tt 1 

ujr r 



2tt 1 
tog T 



r 



If the integral expressions for ui r and the r dependent 
parts of r and Y^ degenerate to elliptic or elementary 



type, i.e. if we consider light or R possesses multiple 
zeros, we can find analytical expressions for ([66|. (|68|l . 
and with the techniques presented in If R has 
only simple zeros and 6=1, the integral ui r is an en- 
try of the fundamental period matrix uj which enters in 
the definition of the period lattice T of the holomorphic 
differentials dz, T = {ojv + co'v \ v, v' £ 1?}. The more 
complicated integrals involving R in (|68p and (f6^|) can be 
rewritten in terms of periods of the differentials of second 
kind dr and of third kind dP(x\,X2) by a decomposition 
in partial fractions. In this way, expressions for u> r , T and 
which are totally analogous to the elliptic case can be 
obtained. 

It follows that the periastron shift is given by 



and the Lense-Thirring effect by 



A 



Lcnsc— Thiiring 



fi<4 — Qg = [Yet, — 



2n 

LOg 



(71) 



(72) 



Another way to access information encoded in the orbits 
is through a frequency decomposision of the whole orbit 
J26(. This will be analyzed elsewhere. 



VII. SUMMARY AND OUTLOOK 

In this paper we derived the analytical solution for 
both timelike and lightlike geodesic motion in Kerr-(anti- 
)de Sitter space-time. The analytical expressions for the 
orbits (r, 8, (f>, t) are given by elliptic Weierstrass and hy- 
perellliptic Klcinian functions in terms of the normalized 
Mino time. We also presented a method for solving dif- 
ferential equations of hyperelliptic type and third kind, 
and applied it to the equations of motion for (f> and t. We 
classified possible types of geodesic motion by an analy- 
sis of the zeros of the polynomials underlying the 9 and 
f motion and discussed the influence of a non- vanishing 
cosmological constant on the orbit types. Some particu- 
lar interesting orbits not present in Kerr space-time were 
shown and a systematic approach for determining the 
last stable spherical and circular orbits was presented. In 
addition, we derived the analytic expressions for observ- 
ables connected with geodesic motion in Kerr-de Sitter 
space-time, namely the deflection angle for escape orbits 
as well as the orbital frequencies, the periastron shift, 
and Lense-Thirring effect of bound orbits. 

The results of this paper can be viewed as the start- 
ing point for the analysis of several features of geodesies 
in Kerr-de Sitter space-time not treated in this publi- 
cation. Although mathematical analogous to the case 
of slow Kerr-de Sitter a complete discussion of orbits in 
fast and extreme Kerr-de Sitter space-times may lead to 
special features and should be carried through. Also, it 
would be interesting to study bound geodesies crossing 
f = (and mayby also the Cauchy horizon for positive 
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f) in general and, in particular, their causal structure. 
In this context the analysis of closed timclike trajecto- 
ries is also of interest. In addition, we not yet consid- 
ered geodesies lying entirely on the axis 6 = 0, 7r or even 
crossing it. Until now, we only considered the Boyer- 
Lindquist form of the Kerr-de Sitter metric which is not 
a good choice for considering geodesies which fall through 
a horizon. Therefore, for a future publication it would be 
interesting to use a coordinate-singularity free version of 
the metric. 

The methods for obtaining analytical solutions of 
geodesic equations presented in this paper are not only 
limited to Kerr-de Sitter space-time. Indeed, they had 
already been used to solve the geodesic equation in 
Schwarzschild-de Sitter [3, EH and Reissner-Nordstrom- 
de Sitter [ItJ space-times. Also the geodesic equations 
in higher-dimensional static spherically symmetric space- 
times 17] were solved by these methods. The same type 



of differential equations is also present in the Plebariski- 
Demiahski space-time without acceleration, which is the 
most general space-time with separable Hamilton- Jacobi 
equation. The analytical solution of the geodesic equa- 
tion in this space-time is given in [ljj but will be elabo- 
rated in an upcoming publication. It will also be inter- 
esting to apply the presented methods to higher dimen- 
sional stationary axially symmetric space-times like the 
Myers-Perry solutions. 

The same structure of equations we solved in this pa- 
per is also present in the geodesic equation of the effective 
one-body formalism of the relativistic two-body problem. 
The effective metric in this formalism can be described 
as a perturbed Schwarzschild or Kerr metric, where the 
pertubation is given in powers of the radial coordinate 
r [27M29| . Therefore, we expect that the polynomial ap- 
pearing in the resulting equations of motion will have 
a higher degree than the corresponding polynomial in 
the Schwarzschild or Kerr case and, thus, that it is nec- 
essary to generalize the elliptic functions used in these 
cases to hyperclliptic functions used in this paper. A 
similar situation can be found in the expressions of ax- 
isymmetric gravitational multipole space-times. For ex- 
ample, some types of geodesies in Erez- Rosen space-time, 
which reduces to the Schwarzschild case if the quadrupole 
moment is neglected, were already solved analytically 
[30l l3l| . We expect that the methods presented in this 
paper will be helpful to solve geodesies in space-times 
with multipolcs. 

Analytic solutions are the starting point for approxi- 
mation methods for the description of real stellar, plan- 
etary, comet, asteroid, or satellite trajectories (see e.g. 
[3^). In particular, it is possible to derive post-Kerr, 
post-Schwarzschild, or post-Newton series expansions of 
analytical solutions. Due to the, in principle, arbitrary 
high accuracy of analytic solutions of the geodesic equa- 
tion they can also serve as test beds for numerical codes 
for the dynamics of binary systems in the extreme stel- 
lar mass ratio case (extreme mass ratio inspirals, EM- 
RIs) and also for the calculation of corresponding gravi- 



tational wave templates. For the case of Kerr space-time 
with vanishing cosmological constant it has already been 
shown that gravitational waves from EMRIs can be com- 
puted more accurately by using analytical solutions than 
by numerical integration [l3j . 

Due to the high precission, the analytical expressions 
for observablcs in Kerr-de Sitter space-time may be used 
for comparisons with observations where the influence of 
the cosmological constant might play a role. This could 
be the case for stars moving around the galactic cen- 
ter black hole or binary systems with extreme mass ra- 
tios where one body serves as test-particle. For example, 
quasar QJ287 (cp. (33|) could be a candidate for observ- 
ing the effects of a non-vanishing cosmological constant. 
In this context it would also be interesting for a future 
publication to derive post-Kerr, post-Schwarzschild, or 
post-Newton expressions for observables. 
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Appendix A: Integration of elliptic integrals of the 
third kind 



In this appendix we will demonstrate the details of the 
integration method of the elliptic integrals of third kind 
which appear in the 9 dependent part of the 4> and t 
motion f|52[) , (|60|) and in the special case of the r motion 
where R has a double or triple zero. We will explain the 
procedure for the example of the integral Ig in ([52]) . As 
this integral is only elliptic if has only simple zeros we 
assume in the following that this is the case. 

Before we demonstrate the calculational steps, we will 
summarize them for convenience: 

1. Cast the expression under the square root in the 
standard Wcierstrass form 4y 3 — giy — gz for some 
constants (72,33, the so-called Weierstrass invari- 
ants. 

2. Decompose the integrand (without the square root) 
in partial fractions. 

3. Substitute y = p(v). 

4. For every partial fraction, rewrite the integrand in 
terms of the Weiertsrass p (double pole) or £ func- 
tion (simple pole). 

5. Integrate the Weierstrass p and Q functions and 
assemble all parts. 
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Let us start from cq. ([54 



/e =2 



a — T> — av' 

a„/(i-i/)vTO; 



:dl/. 



where we assumed that the original integration path was 
fully contained in (0, With the substitutions v = 
and £ = ^- (4y — ^-). as in subsection IIV Al we obtain 



2J So A e (e-l)v^ 



|03| 

2a 3 



((a - P)(V - f ) - aasXV - f )dj/' 



/3 



53 

(Al) 



where is defined in (|40[) , #2 , ff3 are defined as in (|42p , 
A e = £ + a 2 A, and Aj, = d x y+d 2 = ±y-^+a 2 A. Now 
we simplify the integrand in (|A1[) (without the square 
root) by a partial fraction decomposition 



((a-V)(4y-f)-aa 3 )(4y-f) 



a 3 A y (4y 
(a-V) 



(1-2 

3 



03) 



03 

4 X 



a 3 A(x - ahV) 



V 



V - d 3 



(A2) 



where d 3 = f§ + if • With the substitution y = p(v) 
we can get rid of the square root in (|Aip as p'(u) = 
± ^/4p 3 (u) — g 2 p(v) — g 3 where the sign has to be chosen 
according to the sign of p' and the branch of the square 
root. The function p'(v) is negative for v G [0,wi] and 
positive for v G [W2, ^ + w i] where 2cJi G R and 2^2 G (D 
are the fundamental periods of p. As 8 = ~ corresponds 
to y = 00 we will have u G [0, u>x] in most cases. 
Altogether, the integral Ig now reads 



4 X 



a 3 A 



(x - aAV)dv' 



(«0 + 



2W 



p(i/) - d 3 



(A3) 



The second and third integral are of third kind because 

fx(v) = (p(v) + ^ and f 2 (v) = (p(v) - da)" 1 have 
simple poles. We will rewrite now f\ and f 2 in terms the 
Weierstrass ^-function, which has a simple zero in 0. The 
reason is, that £ can easily be integrated in terms of the 
Weierstrass a-function 



((v')dv' = loger(v) — logcr(i> ) 



(A4) 



We only demonstrate the procedure for fx which is to- 
tally analogous to the procedure for f 2 . The function fx 
has two simple poles vx and v 2 lying in the fundamental 
domain {2au>i + 2buj 2 | a, b G [0, 1)}, where u>x and u 2 as 



above, with p(vx) = — 4 



and p(t>) + ^ in neighbourhoods of Vi yields 



p{v 2 )- An expansion of fx 
)ods of Ui yields 

fx(v) — a~x,i{v — Wi) -1 + holomorphic part , 

(A5) 

p(v) + — - = p'(fi)(u — «i) + higher order terms (A6) 
dx 

for some constants a_i : i. Now a comparison of coeffi- 
cients gives 



l = /i(«) + £ 

= a-i,ip'(u») + higher order terms , 
that is a_ M = (p'{vi))~ l . (A7) 

It follows that the function fi(v) — Y2i=x ^p'Jv ) 1 1S an 
elliptic function without poles and, therefore, equal to a 
constant A which can be determined by = /(0). This 



yields A =-Y,Ux^ 

A(») = E 



and 



(=1 



p'OO 



(A8) 



Note that p'{v 2 ) = p'{2ojj - vx) = p'{—vx) = —p'(vx)- 
The expression for p'(«i) can be determined by the differ- 
ential equation p'(x) = ±^ Ap(x) 3 — g 2 p(x) — g 3 , where 
again the sign of the square root has to be chosen accord- 
ing to the sign of p' . An intergation of fx yields 



f fx{v')dv' = j^-±- 



C(vi)(v-v )+log(a(v-Vi)) 



log(cr(w - Vi)) 



(A9) 



In the same way we can integrate f 2 , where we assume 
that i>3 and V4 are the simple poles of f 2 in the funda- 
mental domain with p(v 3 ) = d 3 = p(vx). Summarized, 
Ie is given by (cp. ([55]) ) 



E 



°3 ( H w s . logjajv - Vi)) . „ 

({Vi)(v - Vq) + - — — rr + 2-Kik, 



Jl 4xp'(w 4 ) V log(cr(u - Vi)) 

■ (a 3 A(x - aAV)(8ix + S i2 ) + V{5 l3 + 5a)) | , 

where v = v(-y) = 27 — j$ t - m by (j43|) and vq = ^(70). The 
integers fcj correspond to different branches of log. 



Appendix B: Integration of hyperelliptic integrals of 
the third kind 



In this appendix we will demonstrate the details of the 
integration method of the hyperelliptic integrals of third 
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kind which appear in the r dependent part of the <j) and t 
motion (|52|) . ([60]) . We will explain the procedure for the 
example of the integral I r in (152p . As this integral is only 
hyperelliptic if we consider timelike geodesies, i.e. S — 1 
and R has only simple zeros we assume in the following 
that this is the case. 

Before we demonstrate the solution steps, we will sum- 
marize them for convinience: 

1. Cast the expression under the square root in the 
standard form u 5 + J2t=o CiVjl ^ or somc constants 

a. 

2. Decompose the integrand (without the square root) 
in partial fractions. 

3. If integrals of first or second kind are present, 
rewrite them as functions of 7. 

4. Rewrite the integrals of third kind in terms of the 
canonical integral of third kind dP(xi, 22). 

5. Rewrite the canonical integrals of third kind in 
terms of the Kleinian sigma functions, such reduc- 
ing them to functions of integrals of first kind. 

6. Express the integrals of first kind in terms of 7 and 
assemble all parts. 

Let us start with eq. (|56|) 



a (f' 2 + a 2 - aD) df' 



which can analogously to section HVBI be transformed to 
the standard form by f = ±l/u + with a zero f R of R 
and get 



-a 



u (±±+f R ) 2 + d(d~V)du 
A f=±1/u+fSa \/u- 6 c 5 R u u 2 



'«o ^-»r=±l/u+rR 1 



"0 



(Bl) 



where R u and C5 arc defined as in section IIVBI and 



-^A u = A F= i +fH , i.e. 



A„ = (u 2 {l-Af 2 R )T2f R Au~A)(u 2 (f 2 R +d 2 )±2f R u+l) 

(B2) 



T u 3 - r R u 4 , 



which is a polynomial of degree 4 in «. Note that for 
geodesic motion, the coordinate f is always contained 
in an interval bounded by two adjacent real zeros of 
the polynomial R or by a real zero and infinity. This 
implies that u = ±(f — T R )~ for a real zero f R of R 
does not change sign on the integration path and, there- 
fore, we can neglect the absolute value of u appearing 



in the integrand if we multiply the hole integral with 
sign(it ) = . Consequently 



k 

a 



U [ f R + «(« - V )\ u2 ± 2f R u + 1 3 



/C5 A-uvRH 



11 du . 
(B3) 



The expression for the integrand in (|B3[) can be sim- 
plified by a partial fraction decomposition 



C5|"0| 



4 



Cn 



u du 



/It 
-0 



Uq V^u 

da 



(u — Ui)^/Rj 



, (B4) 



where u.i, 1 < i < 4 denote the zeros of A u and 
Co, C±, 6*2,1 are quite complicated expressions dependent 
on the parameters and the zero f R of R which may be 
calculated by a Computer Algebra System. 

The first two integrals in this expression are of first 
kind and can be solved analogous to section IIVBI (1481) . 
i.e. 



u udu 

Uq V ' Ru 

u du 



I'll 



'Ru 



£5(7 - 70) , 

du 



(B5) 



du 



uq v R u J 00 v Ru 

= -/(V^STO - 7f,in) + fWc5l - 7r,in) 



(B6) 



where again j f ,m = V^57o + X 



00 u' du' 



'R„ 



with 



it 



i (fo — r R ) only depends on the initial values 70 and 
uq, and / describes the (9-divisor, i.e. a ((/(z), z)*) = 0. 

The four integrals in (|B4[) containing (it — it;) -1 are 
in general of third kind and can be expressed in terms 
of the canonical integral of third kind J dP{x\, X2 ) - The 
most simple construction of a differential of third kind 



dP(xi,x 2 ) = 



V + Vi _ V + V2 \ dx 
x — x\ x — X2 ) 2y 



has simple poles in the points (x\, yi) and (^2, 2/2) of the 
Ricmann surface of y 2 = g{x), g a polynomial, with resid- 
ual +1 and —1, respectively (cp. |34l. l35j). In particular, 
we get 



dP(uf,uT), (B8) 




where 



is the pole Ui located on the pos- 



itive branch of the square root and u~ = (ui, — R Ui ) is 
the pole Ui located on the negative branch of the square 
root. Based on Riemann's vanishing theorem (see e.g. 
1341) the canonical differential of third kind can be ex- 
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pressed in terms of Kleinian a functions by 



1 a(f u dz-2[ u * dz) 
dP(u+,ur) = - log ho_^ 

2 a(f u dz-2f u ' dz) 



1 a(f Uo dz~2 f' dz) 

- log Uo ° ^2£_ i . 

2 a(Qdz-2Q dz) V/u„ 



(B9) 



where dz is the vector of canonical differentials of the 
first kind and df the vector of canonical differentials of 
the second kind 



dZi = ,_ , i = l,2, 



'Ru 



5-i 



dri = ^2(k + 1 - i) 



Ck+i+i u k du 
c 5 



(BIO) 

i = 1,2. (Bll) 



Finally, we rewrite (|B9p in terms of the affinc parameter 
7. By (|B5|1 and (|B6|) we can express J^ 1 dz as well as the 
arguments of the a functions f" <iz = f" <iz — f°° dz as 
functions of 7. If we define w = w^j) = y / ci7 — 7f ,in and 
wo = ^(70) the integral 7 r is given by (cp. (|59|) ) 



auQ 



C5\U 



C 1 (w-w )+Co(f(w)-f(w )) 



E 



C 2 ,« 



1 q(W+H) 1 a(W+(w )) 

2 ° S a(W-(w)) 2 ° S a(W-(w )) 



(f(w) - f(w ),w- w ) 



df 



where W + {w) := (f(w),w) 1 — 2 J u< dz and W (w) 
(f(w),wY-2£> dz. 
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